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ABSTRACT 

The morphology of the total gas and the neutral hydrogen (H i) distributions during 
the cosmological epoch of reionization can be quantified with Minkowski Functionals 
(MFs) of isodensity surfaces. We compute the MFs from the output of a high-resolution 
numerical simulation which includes explicit treatment of transfer of UV ionizing ra- 
diation. "Galaxies" identified in the simulation using semi-analytic models of galaxy 
formation are assumed to be the sole sources of UV photons. The MFs of the total 
gas distribution are well described by the analytic expressions derived for lognormal 
random fields. The statistical properties of the diffuse H I depend on the gas distri- 
bution and on the way ionized regions propagate in the inter-galactic medium (IGM). 
The deviations of the MFs of the H I distribution from those of a lognormal random 
field are, therefore, caused by reionization. We use the MFs to discriminate between 
the various stages of reionization in the simulation. We suggest a simple model of 
reionization which reproduces the MFs derived from this simulation. Using random 
realizations of lognormal density fields, we also assess the ability of MFs to distin- 
guish between different reionization scenarios. Our results are relevant to the analysis 
of future cosmological twenty-one centimeter maps. 
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1 INTRODUCTION 

Current observations of the Universe probe many critical 
stages in its history. Snapshots of initial density fluctuations 
similar to those which have led to today's observed struc- 
ture are taken by cosmic microwave background (CMB) ex- 
periments. Galaxy redshift surveys, and QSO spectra probe 
the Universe from the present epoch back to high redshifts 
(z < 6). These data have substantially tightened our grip on 
the cosmological model and its fundamental parameters (e.g. 
Spergel et al. 2003; Percival et al. 2002; Zehavi et al. 2002; 
McDonald et al. 2005; Croft et al. 2002; Nusser & Haehnelt 
2000; Desjacques & Nusser 2005; Zaroubi et al. 2005; Tytler 
et al. 2004). But most data leave out two important con- 
secutive epochs. The first is the Dark Ages starting after 
the formation of atomic hydrogen 380,000 years after the 
Big-Bang and ending with the onset of galaxy and, possi- 
bly, QSO formation. The second is the epoch of reionization 
(EoR) during which the diffuse hydrogen, the dominant pri- 
mordial element in the Universe, is reionized by radiation 
emitted from the emergent luminous objects. The end of 
the EoR, defined when most of the diffuse hydrogen is ion- 



ized, must have occurred before the Universe is a Gyr old, as 
indicated by the amount of QSO light absorbed by neutral 
hydrogen (H I) (e.g. Gunn & Peterson 1965; Becker et al. 
2001). 

The EoR leaves several observational signatures. Thom- 
son scattering of CMB photons off electrons released by 
reionization has two consequences. First, it damps primary 
anisotropies by a factor that depends on the optical depth for 
Thomson scattering. Second, it introduces secondary tem- 
perature and polarization anisotropies resulting from vari- 
ations in the density and velocity fields of free electrons in 
the IGM. Secondary anisotropies, expected to be observed 
by ground based interferometers, will contain information 
on line-of-sight integrated quantities related to the ionized 
gas in the IGM. They are, however, insensitive to the de- 
tails of the H I component and its temporal evolution. A 
probe of these details can be found in QSO spectra. Light 
emanating from QSOs active during the EoR is absorbed by 
the intervening H I. QSO spectra of the most distant QSOs 
are therefore a direct probe of the H I line-of-sight distri- 
bution. Currently, spectra of several QSOs at sufficiently 
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high redshifts (z ~ 6) have been observed (e.g. Becker et al. 
2001). However, bright QSOs are hard to detect at high red- 
shifts, and are therefore mainly suitable for probing the late 
stages of reionization. The inference of constraints on the 
three-dimensional (3D) structure of reionization from QSO 
spectra can be very tricky (e.g. Nusser et al. 2002). Fur- 
ther, even mild H I densities can completely obscure QSO 
light at the relevant frequencies, preventing an accurate es- 
timate of the local ionized fraction. A promising probe of 
reionization, especially of its initial stages, is maps of the 
redshifted 21-cm line. These maps will contain 3D infor- 
mation on the distribution of H I and the ionized fraction 
as a function of time. The 21-cm line is produced in the 
transition between the triplet and singlet sub-levels of the 
hyperfine structure of the ground level of neutral hydrogen 
atoms. This wavelength corresponds to a frequency of 1420 
MHZ and a temperature of T* = 0.068 K. The spin temper- 
ature, T s , is defined according to the relative population of 
the triplet to the singlet sub-levels. An H I region would be 
visible against the CMB either in absorption if T s < T CBR 
or emission if T s > T CBR , where T CBU « 2.73(1 + z) K 
is the CMB temperature. There are various mechanisms 
for raising T s significantly above T CBR during the EoR and 
hence a significant cosmological 21-cm signal is expected. 
The prospects for measuring the 21-cm signal are excellent 
in view of the various radio telescopes designed for this pur- 
pose. Despite foreground contaminations (e.g. Oh & Mack 
2003; Di Matteo, Ciardi & Miniati 2004), the collective data 
obtained by telescopes like the Low Frequency Array (LO- 
FAR, e.g. Kassim et al. 2004) , the Primeval Structure Tele- 
scope (PAST, e.g. Pen, Wu & Peterson 2004; Peterson, Pen 
& Wu 2006), the Square Kilometer Array (SKA, e.g. Carilli 
2004a & 2004b), and the Mileura Widefield Array (MWA, 
Bowman, Morales & Hewitt 2005) should, at the very least, 
help us discriminate among distinct reionization scenarios. 

The statistical properties of the H I distribution are dic- 
tated by the total gas density field and the propagation of 
ionizing radiation (whether UV or X-rays) emitted by galax- 
ies and gas accreting black holes (e.g. Ricotti & Ostriker 
2004). Therefore, 21-cm maps contain a wealth of cosmolog- 
ical information. However, extracting this information can 
be cumbersome. Correlation analysis of 21-cm maps measure 
the power of the signal as a function of scale and may con- 
strain the underlying mass power spectrum. But correlations 
are sensitive neither to non-Gaussian features nor to the 
morphology of the large scale H I distribution. Minkowski 
functionals (hereafter MFs, Hadiwiger 1957), however, pro- 
vide a complete set of measures of the topology of any struc- 
ture. They have been used to study the topological and 
geometrical properties of isodensity surfaces of the galaxy 
distribution (e.g. Mecke et al. 1994; Schmalzing & Buchert 
1997; Kerscher et al. 1997, 1998 & 2001). MFs have also 
been used in the analysis of CMB anisotropics , in particular 
for quantifying non-Gaussian features (Schmalzing & Gorski 
1998; Novikov, Schmalzing & Mukhanov 2000; Schmalzing, 
Takada & Futamase 2000; Shandarin et al. 2002; Komatsu 
E. et al. 2003; Rath & Schuecker 2003; Eriksen et al. 2004). 
In this paper we propose the MFs of H I isodensity surfaces 
as measures of the morphology of reionization. We expect 
the MFs to be most suitable for the analysis of 21-cm maps 
which will provide a direct probe of H I. We focus on the 



ability of MFs to distinguish between the various stages of 
reionization and also between different reionization scenar- 
ios, e.g. by an X-ray background, by UV ionizing patches 
preferentially either in voids or dense regions. We use a high- 
resolution numerical simulation (Ciardi, Ferrara & White 
2003) in which reionization is assumed to be driven by UV 
photons emitted by "galaxies" identified in the simulations 
using a semi-analytical model of galaxy formation. An ex- 
plicit Monte-Carlo scheme has been used to trace the pho- 
tons in the simulation box. The simulation has a relatively 
small box size and it assumes only stellar type sources for 
the ionizing radiation. To overcome these limitations, we use 
an analytical model of reionization which produces a mor- 
phology similar to the one obtained from the simulation. We 
then employ simple reionization recipes (e.g. Benson et al. 
2001) in random realizations of lognormal density fields to 
explore the H I morphology in different reionization scenar- 
ios. Although this may seem simplistic, it allows us to assess 
the ability of MFs to discriminate among these scenarios. 

The paper is organized as follows. In § 2 we review the 
definition of MFs. In § 3 we examine the morphology of 
the gas and of its neutral component in a high resolution 
simulation of the IGM. In § 4 we present MFs from different 
reionization scenarios. We conclude with a summary and a 
discussion of the results in § 5. 



2 THE MINKOWSKI FUNCTIONALS 

Let u{x) denote a scalar function of position, x, defined in a 
3D volume, V. Assume that u has zero mean over the volume 
V and let a be the r.m.s value given by a 2 = (u 2 ) . Given a 
threshold value u thr we consider the excursion set F v formed 
by all points, x, with u(x) > va, where v — u thr /a. The 
global morphology of u(x) is then described by calculating 
the MFs as a function of w thr . In 3D space there are four 
MFs, namely 

V (u) = iy d 3 xQ{uv-u(x)), (1) 
= (2) 

V " {V) = 6^V J 9F ds[Kl{x) + K2{x)] > (3) 
Vz(v) = I ds k-l{x)k2{x), (4) 

J dF u 

where O is the Heaviside step function and, ki and K2 are 
the principal curvatures (inverse of the principal radii) at 
a point x on the surface. The first MF, Vo(z/), is the total 
volume of the regions with u above the threshold u thr = va. 
The remaining three MFs are calculated by surface integra- 
tion over the boundary, dF, of the excursion set, where ds 
denotes the surface element on dF. The MF V\ is a measure 
of the surface area of the boundary, Vi is the mean cur- 
vature over the surface, and Vz is the Euler characteristic, 
X- Compared to a sphere, oblate (pancake-like) surfaces are 
characterized by a large surface area, Vi , and a small mean 
curvature, V2. The Euler characteristic, x, can be expressed 
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X [h"« Mpc] X [h" 1 Mpc] X [h" 1 Mpc] 

Figure 1. Fluctuations in the total gas density (left column), the H I density (middle column), and the ionized fraction X H n = n H n /fi„ 
(right column) in a slice of the simulation box at z = 17.64 (top), z = 15.14 (middle row), and z = 13.68 (bottom). The mean mass 
weighted ionized fractions, /i on , corresponding to these redshifts are indicated to the right of the figure. In the panels showing gas and 
H I fluctuations, the thick solid contour indicates the mean hydrogen density n H = 2.66 X 10 — 'cm -3 , while the thin solid and dotted 
contours, respectively, represent densities above and below n H . The contour spacing in these panels is O.ln^. In the panels showing 
X H . , the thick contour shows the value of /i on and the thin solid and dotted contours, respectively, denote fractions above and below 
/i on . The contour spacing in the panels in the right column is 0.1. 



in terms of the genus, g, as \ — 2(1 — g), where a sphere 
has g = and a torus has g — 1 (Coles, Davies & Pearson 
1996; Gott, Weinberg & Melott 1987). The genus can also 
be viewed as the number of distinct ways an object can be 
sliced without being disconnected into disjoint parts. Hence, 
two separate spheres have g = — 1 since they can be consid- 
ered as an object already made of two disconnected parts. 
Note that, while the Euler characteristic, x> is additive, the 
genus is not so (e.g. Coles, Davies & Peterson 1996). 



For a Gaussian random field in 3D, the average MFs per 
unit volume have the following analytic form (Schmalzing & 
Buchert 1997) 



© 0000 RAS, MNRAS 000, 000-000 



4 Gleser et al. 




Figure 2. MFs of excursions sets defined according to a threshold in the field u = ln(l + <5 gaa ) — (dn(l + <5 gas )y, where 5 ga3 is the gas 
density contrast smoothed with a Gaussian window of a width of 0.5 h _1 Mpc. The squares are obtained from the gas distribution in the 
simulation at rcdshift z = 13.68. The thick dotted lines represent the analytic expressions for Gaussian fields with A estimated from the 
field u. The thin solid lines are the MFs of ten Gaussian random fields realization with the power spectrum of the ACDM cosmology. 



V2{v) = l^ VeXP (-^)' (5C) 
W = -|=(z, 2 -l)exp(-lz/), (5d) 

where A = ^/aJ/W 2 , a = (u 2 ) 1 ' 2 , o x = (|Vu| 2 ) 1/2 . 

Given a sampling of the field, u(x), on a discrete cubic 
lattice we use two methods to calculate the MFs (Schmalz- 
ing & Buchert 1997). The first method is based on Koen- 
derink invariants (Schmalzing 1996; Koenderink 1984; ter 
Haar Romeny et al. 1991). In this method, the local curva- 
tures are expressed in terms of geometric invariants formed 



from the first and second derivatives of the field. The sec- 
ond method is based on Crofton's formula (Crofton 1868; 
Schmalzing & Buchert 1997), which requires counting of the 
number of grid points, edges, surfaces and cells within the 
excursion set. The appendix A contains a more detailed de- 
scription of these methods. We have compared the MFs ob- 
tained from the application of the two methods on the gas 
and H I density fields in the simulations. When the fields 
are interpolated on a 128 3 cubic grid using a simple cloud- 
in-cell algorithm, we find a difference of more than 50% be- 
tween the corresponding MFs. The difference between the 
results obtained with the two methods is reduced to < 10% 
by smoothing the density fields with a Gaussian filter of 
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Figure 3. The MFs of gas (thick lines) and H I (thin lines) isodensity surfaces in the simulation. The gas and H I density fields are 
smoothed with a Gaussian window of 0.5 h _1 Mpc in width. Shown are results for redshifts, z = 17.64 (solid lines), 15.14 (dotted lines), 
and 13.68 (dashed lines) corresponding to mass weighted mean ionized fraction of /; on = 0.08, 0.5, and 0.87, respectively. 



0.5 h _1 Mpc width. Therefore, we work with fields smoothed 
with a filter of 0.5 h _1 Mpc. This filtering roughly corre- 
sponds to the instrumental resolution of LOFAR (Kassim 
et al. 2004). We have examined the difference between the 
methods using a finer grid of 256 3 and obtained similar dif- 
ferences between the methods. The Koendernick invariants 
method, however, seems to be more robust and we only 
present results obtained with this method. 



3 SIMULATIONS OF IGM REIONIZATION 

The simulation we adopt has been run by Ciardi, Ferrara 
& White (2003). This is a high- resolution numerical sim- 
ulation that follows the evolution of the dark matter and 
diffuse gas distribution. It is combined to a semi-analytic 
model of galaxy formation which tracks the sources of ion- 
ization (Kauffmann et al. 1999), and to the Monte-Carlo 
radiative transfer code CRASH which follows the propagation 
of ionizing photons in the IGM (Ciardi et al. 2001; Maselli, 
Ferrara & Ciardi 2003). The simulation also assumes that 
the gas in the IGM consists only of hydrogen. The parame- 
ters are those of a ACDM "concordance" cosmology, with 
n m = 0.3, fi A = 0.7, h = 0.7, Q b = 0.04, n = 1 and 
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as — 0.9 (Spergel et al. 2003). The smallest resolved ha- 
los have a mass M ~ 10 9 M®, and start to form at z ~ 20. 
Each simulation output provides a mock catalogue of galax- 
ies which are described, among others, by their position, 
star formation rate and UV emissivity. For the specific run 
considered here (labeled L20 in the original paper) reioniza- 
tion is completed by z ~ 13, and the associated Thomson 
scattering optical depth is r c 0.161, consistent with the result 
from the first year WMAP data (Spergel et al. 2003). The 
total gas number density, n H , and the neutral gas number 
density, n H t , are provided on a cubic grid of 128 3 . Since we 
have found that MFs computed directly from the density 
fields sampled on a grid without any additional smooth- 
ing are rather noisy, we smooth the gas and H I density 
fields with a Gaussian filter of width of 0.5 h _1 Mpc. All re- 
sults will be presented for the smoothed fields. Using these 
fields, we compute the ionized fraction, X H u = n H n /n H 
(n H jj = n H — n H t ), and the mass weighted mean ionized 
fraction, /i on = (n H X H n ) v / ( n u) V ' Snapshots at redshift 
z = 17.64, z = 15.14, and z = 13.68 will be used. These cor- 
respond to a mean ionized fraction /i 0n = 0.08, /i on = 0.5, 
and /ion = 0.87, respectively. 

Fig. 1 shows contour maps of the smoothed gas and 
H I density fluctuations (left and middle columns, respec- 
tively) at z = 17.64 (top row), z = 15.14 (middle row), 
and z = 13.68 (bottom row). The figure also shows the ion- 
ized fraction, X H u (right column). In the panels to the left, 
the increase in the contour levels as we go from the top 
to bottom is the result of the cosmological growth of fluc- 
tuations (e.g. Peebles 1980). Non-Gaussian features are al- 
ready significant at those redshifts as indicated by the lack of 
symmetry between the distribution of regions with densities 
above and below the mean. A comparison between the left 
and middle columns shows significant differences between 
the morphology of the gas and H I distribution already at 
the initial stages of reionization when the mean ionized frac- 
tion is only /i 0n = 0.08. At these stages, the maps of A H 
(right column) indicate that the ionized regions are mainly 
individual nearly spherical patches. Recombinations play an 
important role at z — 17.64 and 15.14 as inferred from the 
broad distribution of X H u values (see also Ciardi & Madau 
2003). 



3.1 The morphology of the gas distribution 

The evolved gas density field tends to have a lognormal prob- 
ability distribution function (e.g. Theuns et al. 1998). Con- 
sequently, we expect the MFs of the log of the gas density 
field to be well approximated by the analytic expressions for 
Gaussian random fields (eqs. (5a), (5b), (5c) & (5d)). In Fig. 
2, we present the MFs of u = ln(l + <5 gaa ) - (ln(l + <5 gaB )) 
where <5 gaB is the smoothed gas density contrast obtained 
from the simulation output at z = 13.68. There is an excel- 
lent match between the MFs in the simulation and the ana- 
lytic expressions. For high density thresholds, the excursion 
sets are formed by individual density peaks each having an 
Euler characteristic equal to that of a sphere. In that regime, 
V3 is therefore proportional to the number density of peaks 
above the threshold. Note that, in the case of a Gaussian 
field, the symmetric form of V3 reflects the symmetry in the 
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Figure 4. The MFs of the neutral hydrogen density as a function 
of the mass weighted ionized fraction from the simulations. We 
present the results for four density thresholds: = O-OOlnu 

(solid lines), ntj = O.Olnjj (dotted lines), njjj = O-lfij] (dashed 
lines),) and = fin (dot-dashed lines), where Sh = 2.66 X 

1CT 7 cm- 3 . 

statistical properties of negative and positive peaks. Similar 
results were obtained at redshifts z = 17.64 and z = 15.14. 
As an estimate of the cosmic variance in the MFs measured 
from the L20 simulation, we compute the scatter around 
the mean MFs drawn from an ensemble of realizations of 
random fields with lognormal probability distribution func- 
tion (PDF). Namely, Gaussian random realizations of the 
concordance ACDM power spectrum are generated in a pe- 
riodic box of size L — 20 h _1 Mpc. These Gaussian fields 
are mapped onto lognormal fields, which are then normal- 
ized to the r.m.s. value of smooth density fluctuations in the 
simulations. The MFs of 10 such lognormal realizations are 
plotted as the thin lines in Fig. 2. The deviations of these 
MFs from the analytic expressions (thick dotted lines) is the 
result of the finite size of the box. The deviations are most 
pronounced in V3 which depends on a~ 3 (see eq. (5d)) and 
are smallest in Vb which does not depends on a at all (see 
cq. (5a)). 



3.2 The morphology of the H I distribution 

The reionization process in the simulation begins with the 
appearance of individual ionized patches generated by the 
propagation of UV ionization fronts originating from the 
earliest galaxies. During a later stage, when /i on ~ 0.5, the 
geometry of ionized regions transforms into a complex struc- 
ture of tunnels and islands. Towards the end of the EoR, UV 
ionizing photons form a nearly uniform background which 
maintains most of the IGM at a local ionization level de- 
termined by the balance of photo-ionization with recombi- 
nations. In Fig. 3, we plot the MFs of the smoothed gas 
(thick lines) and H I (thin lines) number densities, for three 
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redshifts. The deviations of the MFs of the H I distribution 
from those of the gas are already significant at z = 17.64 
when the ionized fraction is only /i 0n = 0.08. Comparison 
between the thin (H I) and thick (gas) lines in the top-left 
panel shows that the highest density regions are all ionized. 
At 2 = 17.64, all regions with density > 1.42n H are ion- 
ized together with some moderate density regions, whereas 
at z — 13.68, all the regions with mean densities and above 
are ionized. Between z = 17.64 and z = 15.14 significant 
reionization occurs in regions of densities below the mean 
(compare solid and dotted thin lines). The MFs at the three 
redshifts are substantially different, making the MFs a good 
discriminator of the various stages of reionization. In Fig. 
4, we examine the evolution of the MFs of the H I distri- 
bution as a function of the mass weighted ionized fraction, 
/ion, for the density thresholds: n\ i \ = O.OOlnH (solid lines), 
h'h i = 0.01n H (dotted lines), n%\ = 0.1n H (dashed lines), 
and n^j = fin (dot-dashed lines), where the mean gas den- 
sity is tih = 2.66 x 10~ 7 cm -3 . As expected, the volume 
functional, Vo{v), decreases with /i on . At the initial stages 
of reionization, the ionized volume of space is made of non- 
overlapping patches, so that the surface area, Vi, is an in- 
creasing function of /ion, and the integrated mean curvature 
functional, Vziy), is negative . When the ionized fraction /i on 
becomes significant, the behaviour of Vi, V2, and the Euler 
characteristic V3, is reversed compared to their behaviour at 
low /ion- The reason is that both the ionized volume in the 
early phase of reionization, and the neutral volume in the 
late stages, are predominantly made of individual patches 
and, therefore, share similar topological properties. 



4 SEMI-ANALYTIC MODELING OF THE 
MORPHOLOGY OF REIONIZATION 

4.1 Modeling the morphology in the simulation 

Our goal here is twofold. First, we will attempt to repro- 
duce the MFs of the H I distribution in the simulation using 
simple recipes that do not require explicit treatment of ra- 
diative transfer. Second, given these recipes we will assess 
the extent to which MFs can discriminate between reion- 
ization scenarios other than those adopted in the particular 
simulation used here. In the recipes A, B, C, and D, that 
we describe below, the H I distribution is obtained by over- 
laying the gas distribution with spherical bubbles devoid of 
neutral gas, while the surrounding gas is kept neutral. In the 
recipes E and F the IGM in the simulation box is ionized cell 
by cell as explained below. All the models assume that a grid 
cell is either completely ionized or neutral. Clearly, this is an 
oversimplification especially in the early stages of reioniza- 
tion (see column to the right in Fig. 1). Yet we will see that 
one of these recipes indeed reproduces the MFs measured in 
the simulations. We now summarize the different recipes. 

• A: spherical bubbles are placed randomly in the simu- 
lation. All bubbles have the same radius, R, independent of 
local gas density. 

• B: the same as A, but with bubble radius, r, depending 
on the local density. We have tried several functional forms, 
and found that the following parametrization gives the best 



match to the MFs of the simulation, 

rM _( R + 6R[l-(l + <5(*)r 1/3 ] 1 + S > (6/7) 3 
[ j I l + <S<(6/7) 3 (6) 

where S(x) is the smoothed density contrast at the position 
x, and R is a constant. 

• C: bubbles of fixed radius R, are added one by one, 
where a newly added bubble is centered on the highest den- 
sity grid point not encompassed by previous bubbles. 

• D: the same as C but with new bubbles centered at 
lowest density point. 

• E: grid cells are ionized according to their density, start- 
ing from the cell with the highest gas density. 

• F: the same as E but starting from the lowest density 
cell. 

We apply each of the recipes described above to the gas 
distribution in the L20 simulation. To facilitate the compar- 
ison with the MFs derived from the H I distribution in the 
simulation, the number of bubbles in A-D, and the number 
of ionized cells in E & F are tuned so that the mean mass 
weighted ionized fraction matches the value measured in the 
simulation at a given redshift. In Fig. 5, we compare the MFs 
obtained with these recipes to those of the actual H I dis- 
tribution in the simulation. Left, middle, and right columns 
show, respectively, results for z = 17.64, 15.14, and 13.88. 
In each of the models A-D, the parameter R is tuned to 
the value giving the best match to the curve of Vo obtained 
from the simulation. We have found that this value mainly 
depends on redshift: R = 0.94, 1.56, and 2.03 h _1 Mpc co- 
moving for z — 17.64, 15.14, and 13.88, respectively. The 
figure shows results for these values of R only. Recipe B 
(thick dashed lines) yields the best fit to the MFs of the 
actual H I distribution in simulation (think solid lines). In 
this model, the radius of bubbles in regions with overden- 
sities < 0.63 is practically zero. Therefore, these regions 
are only ionized by sources in nearby denser regions. This 
reflects the relative lack of ionizing sources (galaxies) in the 
voids of the simulation volume. 

The numerical simulation we use has a relatively small 
box (20 h _1 Mpc) so that cosmic variance may introduce 
large uncertainties in the estimates of the MFs (see thin lines 
in Fig. 2) . Furthermore, it is assumed in this simulation that 
reionization is caused solely by UV radiation. The amount of 
UV photons emitted by a "galaxy" in the simulation is sen- 
sitive to the adopted parameters of the semi-analytic galaxy 
formation model. The distribution of the "galaxies" , which 
greatly affects the development of the early stages, also de- 
pends on the galaxy formation model. Exploring the param- 
eter space of the galaxy formation model and all possible sce- 
narios with large volume simulations is currently unfeasible. 
Nonetheless, Figs. 5 demonstrate that the MFs can distin- 
guish between the different reionization scenarios examined 
here. 

To reduce the effect of cosmic variance in the MFs we 
apply the recipes A-F to lognormal random fields in a cu- 
bic box of 50 h _1 Mpc on the side (more than 15 times the 
volume of the L20 simulation). The corresponding MFs are 
shown in Fig. 6. The larger box size significantly reduces the 
noise in the MFs, especially for V3. In general, there are only 
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-1.0 -0.5 0.0 0.5 -0.5 0.0 0.5 -0.5 0.0 0.5 

lo g( n thr/ n H ) lo g( n thr/n H ) lo g( n thr/ n H ) 

Figure 5. Comparison of the MFs of the smoothed H I density field from the simulation with those obtained from the recipes A, B, 
C, D, E, & F. Left, middle and right columns correspond to z = 17.64, z = 15.14, and z = 13.68, respectively. The thick solid lines 
represent the MFs from the H I distribution in the simulation. The thick dotted corresponds to recipe A, thick dashed to B, thin solid to 
C, thin dashed to D, thin long-dashed to E, and dash-dotted to F. The best match to the simulation is obtained with B (thick dashed). 
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lo g( n thr/ n H ) lo g( n thr/n H ) lo g( n thr/ n H ) 

Figure 6. The same as the previous figure, but with recipes A-F applied to a lognormal random density field in a cubic box of 50 h -1 Mpc 
on the side. Here also, recipe B (thick dashed) yields the best match to the MFs from H I field in the simulation (thick solid). 
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Figure 7. MFs of excursions sets obtained from the H I distri- 
bution in an X-ray pre-reionization scenarios. Plotted arc curves 
for three values of X u at mean gas density, as indicated in the 
figure. All curves are for z = 17.64. 



Figure 8. The MFs from the previous figure plotted against 
threshold in v = [ln(n H : ) — ^ln(n H /cr where cr 2 = 

^ [ln(n R j) — ^ln(n H The analytic MFs for Gaussian ran- 

dom field arc shown as the dash-dotted lines. 



minor differences between Figs. 5 and 6, indicating that the 
lognormal approximation yields results in agreement with 
the simulation. This is expected since the morphology of 
the total gas density is well approximated by the MFs of 
a lognormal field. Therefore, the morphology of reionization 
can also be modeled using lognormal fields without resorting 
to detailed simulations. 



4.2 Morphology in an X-ray pre-reionization 
scenario 

Reionization by UV sources can be preceded by a pre- 
reionization epoch in which the gas is partially ionized by 
an X-ray background produced by accretion onto mini black 
holes (Ricotti & Ostriker 2004). The mean free path of X-ray 
photons is much larger than both the correlation length and 
the mean separation of the sources. Therefore, in this sce- 
nario, X-ray photons act as a uniform ionizing background 
and the H I fraction, X u 1 = n H 1 /n H , at a point with den- 
sity contrast 5, is determined by the equation for photo- 
ionization equilibrium, 

rX HI =4 2) I fi(l + ( 5)(l-X HI ) 2 , (7) 
where F is the X-ray photo-ionization rate per hydrogen 

(2) 

atom, qJj j is the recombination cross section to the second 
excited atomic level, n is the mean density of total (ion- 
ized plus neutral) hydrogen, and X H is the H I fraction. 
According to Ricotti & Ostriker (2004), the mean H I frac- 
tion should be maintained at a level X H : ~ 0.8 — 0.9 in 
the redshift range 20 < z < 40 in order to yield a large op- 
tical depth t ~ 0.1. This level of reionization is ideal for 
producing a significant 21-cm signal (Nusser 2005a; Kuhlen, 
Madau & Montgomery 2006). In order to model the MFs 
of H I isodensity surfaces in this scenario, we use the log- 



normal gas density fields which have been computed from 
Gaussian random fields in a box of 50 h _1 Mpc on the side. 
We fix the ratio a H \n/V by assuming that X H I is given at 
mean density, i.e. at 8 = 0. Then, for any 8, X H : can be 
easily found using eq. (7). In Fig. 7 we present the MFs of 
excursions sets obtained from X-ray background H I fields 
at redshift z = 17.64 for values of the H I fraction at mean 
density of X u 1 = 0.95 (solid lines), 0.9 (dotted lines), and 
0.8 (dashed lines). Changing the H I fraction shifts the MFs 
horizontally, but does not affect their overall shape. This 
is not surprising since all these models have similar topolo- 
gies. In Fig. 8, we plot the MFs as functions of excursion 
sets obtained from the field ln(n Hl ). Also plotted, as the 
dash-dotted lines, are the analytic expressions for Gaussian 
random field. The figure demonstrates that for the low levels 
of ionization obtained in X-ray pre-reionization, the MFs of 
the H I distribution are close to those of a lognormal field 
which also describes the morphology of the gas distribution 
in our case. 



5 DISCUSSION 

Among all possible measures of global morphology, only 
Minkowski Functionals (MFs) are invariant under rotational 
and translational transformations (Hadwiger 1957). But can 
the MFs distinguish between various stages of EoR or be- 
tween different reionization scenarios? We have assessed this 
question using a high-resolution numerical simulation in 
which reionization is driven by UV photons emitted from 
stellar type sources and by employing simple recipes to de- 
scribe different reionization scenarios. We have found that 
the MFs of isodensity surfaces of the gas distribution in the 
simulation are well described by those of a field having a 
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lognormal distribution function. In a UV rcionization sce- 
nario, the deviations from the lognormal MFs are sensitive 
to the distribution of ionized regions. The MFs of the H I 
distribution depend on a multitude of factors among which 
are: the gas density field and its correlation with the dis- 
tribution of ionizing sources, the abundance of sources, and 
the intensity and type (UV versus x-rays) of the emitted ra- 
diation. Reionization by UV begins with individual patches 
appearing around sources associated with the most massive 
and rarest haloes (e.g. Benson et al. 2005). The topology 
of reionization in these early stages is substantially differ- 
ent from that of later stages when ionized regions form a 
complex network of tunnels in addition to isolated patched. 
This difference in topology is reflected in the shapes of the 
MFs. X-rays have a substantially larger mean free path than 
UV photons. Therefore, x-rays generated by an early gener- 
ation of black holes that tend to form a uniform background 
that may partially ionize the IGM. We have considered the 
MFs of the H I distribution in an epoch of pre-reionization 
caused by an x-ray background. Using a semi-analytic mod- 
eling of this epoch, we find similar H I and gas topologies. 
This is due to the low local ionized fractions obtained in 
these models. 

Despite the complexity of the reionization process as 
given by the simulation (cf. Ciardi, Ferrara & White 2003), 
its topology can be described using simple semi-analytic 
recipes that do not involve a detailed treatment of radia- 
tive transfer. Application of some of these recipes on lognor- 
mal random fields yield a good match to the MFs obtained 
from the simulation. However, one may envisage a variety of 
alternative UV reionization scenarios (e.g. Miralda-Escude, 
Haehnelt & Rees 2000), depending on the details of galaxy 
formation and the dark matter distribution. It therefore re- 
mains to be confirmed whether semi-analytic recipes can 
actually match the MFs obtained from detailed simulations 
of alternative scenarios corresponding to physical assump- 
tions different from those adopted in the simulation we have 
used in this paper. 

In the near future, the next generation of radio tele- 
scopes will provide maps of twenty-one centimeter (21-cm) 
line from neutral IGM and probe the three dimensional dis- 
tribution of H I during the EoR. It will be then possible to 
use MFs to study the morphology of the reionization pro- 
cess as described in this work. However, the observed 21- 
cm maps will be affected by redshift distortions (e.g. Ali, 
Bharadawaj & Pandey 2005; Barkana & Loeb 2005; Nusser 
2005b), instrumental noise, and foreground contamination. 
The usefulness of MFs derived from 21-cm maps which in- 
clude all these effects has yet to be demonstrated. This will 
be the subject of future work in which the modified method- 
ology of Benson et al. (2001, 2005) for studying reionization 
will be implemented in a high resolution N-body simulation 
of a large box (~ 140 h _1 Mpc on the side). 

The background density parameters and the value of as 
of the simulation used in this paper are consistent with the 
first year WMAP data (Spergel et al. 2003). Although the 
third year WMAP data give similar values for the density 
parameters as the first year data, it yields er 8 = 0.74+g'gg 
and t c = 0.09t£'o3 (Spergel et al. 2006). The difference in 
the values for as and r e can amount to substantial changes 



in the reionization history (Benson et al. 2005). This work 
has been completed before the the third year data has been 
published and therefore simulations have been run with cos- 
mological parameters inferred from the first year WMAP 
data. Nevertheless, our aim at this stage is to illustrate the 
use of MFs for studying the history of cosmic reionization. 
Implications of MFs of 21-cm maps generated from a suite 
of simulations will be presented elsewhere. 

Auto-correlations estimated from 21-cm maps have also 
been proposed as a statistical measure of H I distribution. 
The advantage of auto-correlations is that redshift distor- 
tions, noise, and foreground contamination can probably be 
easily incorporated in or removed from them, while these 
effects might play a more important or ambiguous role in 
MFs of 21-cm maps. Nevertheless, as auto-correlations are 
insensitive to morphology (e.g. a Gaussian random field and 
another highly non-Gaussian one, can have identical auto- 
correlation functions), MFs are better discriminators be- 
tween different reionization scenarios. Therefore, an explo- 
ration of MFs of 21-cm maps is worthwhile. 



ACKNOWLEDGMENTS 

AN wishes to thank the Max-Planck-Institut fur Astro- 
physik in Garching for its generosity and support. AN & 
LG acknowledge a useful discussion with T. Buchert. VD 
acknowledges the support of a Golda Meir fellowship. 



REFERENCES 

Ali Sk. S., Bharadwaj S., Pandey B., 2005, MNRAS, 363, 251 
Barkana R., Locb A., 2005, ApJ, 626, 1 
Becker R. H. et al. , 2001, AJ, 122, 2850 

Benson A. J., Nusser A., Sugiyama N., Lacey, C.G., 2001, MN- 
RAS, 320, 153 

Benson, A. J., Sugiyama N., Nusser A., Lacey C. G., 2006, MN- 
RAS, 369, 1055 

Bowman J. D., Morales M. F., Hewitt J. N., 2006, ApJ, 638, 20 
Carilli C. L., Furlanetto, S., Briggs F., Jarvis M., Rawlings S., 

Falcke H., 2004a, NewAR, 48, 1029 
Carilli C. L., Gnedin N., Furlanetto S., Owen F., 2004b, NewAR, 

48, 1053 

Ciardi B., Ferrara A., Marri S., Raimondo G., 2001, MNRAS, 
324, 381 

Ciardi B., Madau P., 2003, ApJ, 596, 1 

Ciardi B., Ferrara A., White S. D. M., 2003, MNRAS, 344, L7 
Coles P., Davics A. G., Pearson R. C, 1996, MNRAS, 281, 1375 
Croft R. A. C, Weinberg D. H., Bolte M., Buries S., Hernquist 

L., Katz N., Kirkman D., Tytlcr D., 2002, ApJ, 581, 20 
Crofton M. W., 1868, Philos. Trans. R. Soc. London, A, 158, 181 
Desjacqucs V., Nusser A., 2005, MNRAS, 361, 1257 
Di Matteo T., Ciardi B., Miniati F., 2004, MNRAS, 355, 1053 
Eriksen H. K., Novikov D. I., Lilje P. B., Banday A. J., Grski K. 

M., 2004, ApJ, 612, 64 
Gott J. R. Ill, Weinberg D.H., Melott A. L., 1987, ApJ, 319, 1 
Gunn J.E., Peterson B.A., 1965, ApJ, 142, 1633 
Hadwigcr H., 1957, Vorlesungen iiber Inhalt, Oberflache und 

Isoperimeterie (Berlin: Springer) 
Kassim N. E., Lazio T. J. W., Ray P. S., Crane P. C, Hicks B. 



© 0000 RAS, MNRAS 000, 000-000 



12 Gleser et al. 



C, Stewart K. P., Cohen A. S., Lane, W. M., 2004, P&SS, 
52, 1343 

Kauffmann G., Colberg J. M., Diaferio A., White S. D. M., 1999, 
MNRAS, 307, 529 

Kerscher M., Schmalzing J., Retzlaff J., Borgani S., Buchcrt 
T.,Gottl6ber S., Muller V., Plionis M., Wagner H., 1997, MN- 
RAS, 284, 73 

Kerscher M., Schmalzing J., Buchcrt T., Wagner H., 1998, A&A, 
333, 1 

Kerscher M., Mecke K., Schmalzing J., Beisbart C, Buchert T., 

Wagner H., 2001, A&A, 373, 1 
Koenderink J. J., 1984, Biol. Cybcrn., 50, 363 
Komatsu E. ct al. , 2003, ApJS, 148, 119 
Kuhlen M., Madau P., Montgomery R., 2006, ApJ, 637, 1 
Maselli A., Fcrrara A., Ciardi B., 2003, MNRAS, 345, 379 
McDonald P. et al. , 2005, ApJ, 635, 761 
Mecke K. R., Buchcrt T., Wagner H., 1994, A&A, 288, 697 
Miralda-Escude J., Haehnelt M., Rees M. J. 2000, ApJ, 530, 1 
Novikov, D., Schmalzing, J., Mukhanov, V. F., 2000, A&A, 364, 

17 

Nusser A., Haehnelt M., 2000, MNRAS, 313, 364 
Nusser A., Benson A. J., Sugiyama N., Laccy C, 2002, ApJ, 580, 
93 

Nusser A., 2005a, MNRAS, 359, 183 

Nusser A., 2005b, MNRAS, 364, 743 

Oh S. P., Mack K. J., 2003, MNRAS, 346, 8710 

Peebles P. J. E., 1980, The Large-Scale Structure of the Universe, 
Princeton University Press 

Pen U. L., Wu X. P., Peterson J., 2004, Chin. J. Astron. Astro- 
phys., submitted (astro-ph/0404083) 

Percival, W. J. et al. , 2002, MNRAS, 337, 1068 

Peterson J., Pen U. L., Wu X. P., 2006, in Kassim N., Perez M., 
Junor M., Henning P., eds, Astron. Soc. Pac. Conf. Scr. Vol 
345, From Clark Lake to the Long Wavelength Array: Bill 
Erickson's Radio Science, Astron. Soc. Pac, San Francisco, 
p. 441 

Rath C, Schuecker P., 2003, MNRAS, 344, 115 

Ricotti, M., Ostriker, J. P., 2004, MNRAS, 352, 547 

Schmalzing J., 1996, Diplomarbcit, 

Ludwig-Maximilian-Universitat, Miinchcn 
Schmalzing J., Buchert T., 1997, ApJ, 482, LI 
Schmalzing J., Gorski K. M., 1998, MNRAS, 297, 355 
Schmalzing J., Takada M., Futamase T., 2000, ApJ, 544, 83 
Shandarin S. F., Feldman H. A., Xu Y., Tegmark M., 2002, ApJS, 

141, 1 

Spergel D. N. et al. , 2003, ApJS, 148, 175 

Spergel D. N. et al. , 2006, ApJ, submitted (astro-ph/0603449) 
ter Haar Romeny B. M., Florack L. M. J., Koenderink J. J., 
Viergever M. A., 1991, in Colchester A. C. F., Hawkc D. J., 
eds, Lecture Notes of Computer Science, Vol. 511, Scale Space: 
Its Natural Operators and Differential Invariants. Springcr- 
Verlag, Berlin, p. 239 
Theuns T., Leonard A., Schaye J., Efstathiou G., 1998, MNRAS, 
303, 58 

Tytler T. et al. , 2004, ApJ, 617, 1 

Zaroubi S., Viel M., Nusser A., Haehnelt M., Kim T. -S., 2006, 

MNRAS, 369, 734 
Zehavi I., 2002, ApJ, 571, 172 



tice grid. Using the techniques pioneered by Koenderink 
(Koenderink 1984; ter Haar Romeny et al. 1991) in two- 
dimensions, Schmalzing (1996) expressed the local curva- 
tures in term of geometric invariants formed from first and 
second derivatives. These invariants are known as the Koen- 
derink invariants. 



2(ii, P u,p) 2 



(Ala) 



(Alb) 



where ki and «2 are the principal curvatures, w(x) is a three- 
dimensional density field in a box of a volume V, u t i and u,ij 
are the first and second spatial derivatives of the field u(x) 
respectively in the i and j directions, and is the Levi- 
Civita tensor. 

Since one can replace the surface integration in eqs. (2), 
(3) & (4) with a spatial mean over the whole volume, the 
MFs become 



v ^ = w 



d xS(i/G — u(x)) \X7u(x) 



(A2a) 



V 2 {v) = --^-— I <i 3 x5{i>cr -u{x))\Vu(x)\[ki{x)+k 2 (x)\, 
6nV Jv (A2b) 



1 f , 

V 3 (v) = — — / d xS(va — u(x))\\7u(x)\ ki(x)k2(x), 

4%V Jv (A2c) 



where v is the density threshold divided by the field r.m.s. 
value, a, and 5 is a delta function. 

The second method is based on Crofton's formula 
(Crofton 1868) which requires counting of the number of 
grid points, edges, surfaces and cells within the excursion 
set. The MFs in three-dimensions are given by 



V (K) = jN s (K), 



V^K) = 



2 1 
9 al 



2 1 



-3N 3 (K)+N 2 (K)], 



V*(K) = [3N 3 (K) - 2N 2 (K) + N^K)] , 



(A3a) 



(A3b) 



(A3c) 



V 3 (K) = — [-N 3 {K)+N 2 {K) - N^K) + N {K)] 



(A3d) 



where the body K is the excursion set of a homogeneous and 
isotropic random field sampled at L points of a cubic lattice 
of spacing a. N 3 (K) gives the number of cubes within the 
body, N 2 (K) is the number of faces and N\(K) is the num- 
ber of edges of the lattice cubes which belong to the body, 
while No(K) gives the number of lattice points contained in 
K. 



APPENDIX A: NUMERICAL CALCULATION 
OF THE MINKOWSKI FUNCTION ALS 

Schmalzing & Buchert (1997) used an integral geometry ap- 
proach and developed two numerical methods to investigate 
the morphology of a sampled density field on a cubic lat- 



in http://physics.technion.ac.il/~lirong/mfs.tar.gz one 
can find the MFs software package for MFs calculation. This 
package includes the Koenderink invariants method, the 
Crofton's formula method and the analytical calculation for 
Gaussian random field. 
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